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Abstract. We report on a previously undocumented range er- 
ror in NASA’s Ice, Cloud and land Elevation Satellite (ICE- 
Sat) that degrades elevation precision and introduces a small 
but significant elevation trend over the ICESat mission pe- 
riod. This range error (the Gaussian-Centroid or “G-C” off- 
set) varies on a shot-to-shot basis and exhibits increasing 
scatter when laser transmit energies fall below 20 mJ. Al- 
though the G-C offset is uncorrelated over periods < 1 day, 
it evolves over the life of each of ICESat’s three lasers in 
a series of ramps and jumps that give rise to spurious ele- 
vation trends of —0.92 to — 1.90cmyr“^, depending on the 
time period considered. Using ICESat data over the Ross and 
Filchner-Ronne ice shelves we show that (1) the G-C offset 
introduces significant biases in ice-shelf mass balance esti- 
mates, and (2) the mass balance bias can vary between re- 
gions because of different temporal samplings of ICESat. We 
can reproduce the effect of the G-C offset over these two ice 
shelves by fitting trends to sample-weighted mean G-C off- 
sets for each campaign, suggesting that it may not be nec- 
essary to fully repeat earlier ICESat studies to determine the 
impact of the G-C offset on ice-sheet mass balance estimates. 


1 Introduction 

NASA’s Ice, Cloud, and land Elevation Satellite (ICESat) 
(Schutz et al., 2005) was an Earth-orbiting laser altimeter 
mission that operated from 2003-2009. ICESat’s primary 
task was to repeatedly measure surface elevations along fixed 
ground tracks over Earth’s polar regions to help quantify 
the contribution of the ice sheets to contemporary sea level 


change. Many studies have used ICESat elevation data to 
estimate volume/mass changes of glaciers (Gardner et al., 
2013), ice shelves (Pritchard et al., 2012), and ice sheets 
(Shepherd et al., 2012), and ICESat data have been com- 
bined with other measurements to increase the spatiotem- 
poral coverage and resolution of surface change estimates. 
These complementary data include airborne laser altime- 
try from NASA’s Operation IceBridge mission (Koenig et 
al., 2010; Kwok et al., 2012; Schenk and Csathd, 2012), 
gravity from the NASA/DLR GRACE mission (Riva et al., 
2009), and elevations from ESA’s ERS-1, ERS-2 and Envisat 
radar altimeters (Zwally et al., 2011; Eiurkmans et al., 2012). 
ICESat will provide benchmark elevations for the planned 
ICESat-2 mission (Abdalati et al., 2010), which is planned 
for launch in 2017 and will extend the satellite laser altime- 
ter record to 17 yr or more. 

Since the ice sheets are so vast, a 1 cm ice-equivalent ele- 
vation change over all grounded ice corresponds to a mass 
change of 134 Gt and a sea level equivalent of 0.37 mm. 
Centimeter-level systematic errors in satellite altimeter mea- 
surements are therefore crucial to ice mass balance esti- 
mates, and ICESat underwent rigorous calibration and val- 
idation to ensure that it would meet its accuracy target of 
1.5cmyr“^ elevation change averaged over 100 x 100 km 
regions (Zwally et al., 2002). ICESat’s post-launch calibra- 
tion of various geolocation parameters was based primar- 
ily on the minimization of elevation residuals from regu- 
larly repeated pointing maneuvers over the open ocean (e.g., 
Luthcke et al., 2005), with additional efforts to calibrate tim- 
ing and geolocation using ground truth (Magruder et al., 
2007, 2010) and to mitigate the impact of detector saturation 
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(Fricker et al., 2005). ICESat post-launch elevation valida- 
tion included crossover analysis to determine ICESat ’s initial 
precision and accuracy (Shuman et al., 2006; Brenner et al., 
2007), followed by long-term elevation comparisons with re- 
spect to stable and/or independently characterized reference 
surfaces (Fricker et al., 2005, Urban and Schutz, 2005; Borsa 
et al., 2007, 2008; Shuman et al., 2009). 

Although ICESat was originally intended to be oper- 
ated continuously throughout its mission (Abshire et al., 
2003), concerns about laser reliability after the failure of the 
first ICESat laser led to it being operated campaign-style, 
whereby data were acquired in a series of 3 3 -day cam- 
paigns spaced 4-6 months apart (Table 1). The ICESat val- 
idation effort focused largely on documenting changes in 
ICESat elevation accuracy from campaign to campaign and 
between different releases of ICESat data (e.g., Fricker et 
al., 2005). Despite ongoing refinements in elevation retrieval 
(reflected in higher product release numbers), multiple stud- 
ies have documented persistent instrument-related elevation 
biases between campaigns (Gunter et al., 2009; Riva et al., 
2009; Siegfried et al., 201 1). More importantly, these “inter- 
campaign biases” exhibit statistically significant (albeit dif- 
ferent) trends over the ICESat mission period (Urban et al., 

2012) . Furthermore, in the absence of a definitive set of in- 
tercampaign biases, researchers estimating ice volume/mass 
balance using ICESat data have taken different approaches 
with respect to intercampaign bias correction, with some 
making no correction (Pritchard et al., 2009; Gardner et al., 

2013) and others applying biases from one of several sources 
(e.g., Gunter et al., 2009; Riva et al., 2009; Zwally et al., 
2011; Shepherd et al., 2012). 

This paper describes a previously unrecognized compo- 
nent of the ICESat intercampaign biases, an inadvertent 
range error (called the Gaussian-Centroid or “G-C” offset) 
that was introduced during the processing of Level 1 data. 
Correcting for this error improves the precision of individual 
elevation measurements and removes a small but significant 
anomalous elevation trend from ICESat data. Using global 
statistics for the GC offset and case studies over the salar de 
Uyuni in Bolivia and two Antarctic ice shelves, we demon- 
strate the potential impact that the correction has on ICESat 
elevation accuracy and ice sheet mass balance estimates. 


2 Data and analysis 

2.1 ICESat campaigns 

Data collection during the ICESat mission took place dur- 
ing 18 separate campaigns between February 2003 and Oc- 
tober 2009 (Table 1). In this paper, we refer to these cam- 
paigns using the standard convention of pairing the number 
of the operational laser with a letter designating each consec- 
utive campaign for that laser (e.g., L2a is the first campaign 
for Laser 2, L3b is the second campaign for Laser 3, etc.). 


Laser 1 operated for only 56 days before it failed and was 
flown only in ICESat’s 8-day exact-repeat calibration orbit. 
Most published studies use data from L2a onwards (after the 
spacecraft had transitioned to its 91 -day repeat orbit) so we 
will focus primarily on Lasers 2 and 3. To clarify the time 
sequence of the laser campaigns, we note that Laser 2 was 
switched off after L2c and then back on again after Laser 
3 failed, which is why campaigns L2d-L2f took place after 
L3k (the final Laser 3 campaign). 

2.2 ICESat data products 

ICESat data are publically available from the National Snow 
and Ice Data Center (NSIDC, http://nsidc.org/data/icesat/ 
index.html). Most users of ICESat elevation data choose one 
of several Level 2 data products containing the geolocated 
positions of individual laser footprints: GLA06 (Global El- 
evation Data), GLA12 (Antarctic and Greenland Ice Sheet 
Altimetry data), GLA13 (Sea Ice Altimetry Data), GLA14 
(Global Land Surface Altimetry Data), and GLA15 (Ocean 
Surface Altimetry Data). Different assumptions about sur- 
face characteristics and the needs of investigators resulted in 
differentiation between these data products. Relevant to this 
study, elevations for GLA06, GLA12, GLA13 and GLA15 
were calculated assuming laser reflection from smooth sur- 
faces such as oceans and ice sheets, which generate simple 
return waveforms that can be represented accurately by a sin- 
gle Gaussian peak. Elevations for GLA14 were calculated as- 
suming complex return waveforms that are best summarized 
by the centroid of the entire waveform, a distinction that will 
be discussed later. 

In our work, we also used the Level lA GLAOl (Global 
Altimetry Data) product to access the individual transmit and 
return waveforms for each laser shot; the Level IB GLA05 
(Global Waveform-based Range Correction Data) product 
for waveform parameters such as centroid, Gaussian fit, and 
skewness; and the Level IB GLA06 product for instrument 
pointing, laser energy, surface reflectance, and various cor- 
rections. Although most users need only Level 2 data (pro- 
cessed geophysical variables) in their work, access to data 
in Level IB (processed instrument) and Level lA (raw in- 
strument) products allows investigators such as ourselves to 
contribute to ICESat calibration and validation efforts. 

2.3 Elevation validation at the salar de Uyuni 

This study arose from a longstanding ICESat validation ef- 
fort using a 45-by-54km reference surface on the salar de 
Uyuni in Bolivia that we first surveyed with kinematic GPS 
in 2002 (Fig. 1) (Fricker et al., 2005; Borsa et al., 2007). Our 
digital elevation model (DEM) of the salar surface showed 
total topographic relief of only 80 cm over 50 km, making 
the salar de Uyuni one of the flattest large natural surfaces 
on Earth. ICESat Track 360 (ascending) and Track 85 (de- 
scending) cross the Uyuni DEM, and over 300 individual 
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Table 1. ICESat campaign metadata and G-C offset statistics. Campaigns are listed sequentially in time and are named as described in the 
text. Laser 2 campaigns are shown in grey to highlight the switching that occurs between lasers during the mission. The valid returns column 
gives the percentage of shots for which a surface elevation was recorded, whieh tends to drop as laser energy deereases. 


Campaign 

Start date 

End date 

# Days 

% Valid 
returns 

# Valid 
returns 

G-C mean 
(cm) 

G-C(t 

(cm) 

Lla/Llb 

2/20/03 

3/29/03 

37 

73% 

85 862 366 

1.77 

1.99 

L2a 

10/13/03 

11/19/03 

37 

63% 

80 770 785 

5.90 

1.83 

L2b 

2/17/04 

3/21/04 

33 

64% 

72 697 059 

0.60 

2.07 

L2c 

5/18/04 

6/21/04 

34 

53% 

62 361334 

-0.16 

3.70 

L3a 

10/3/04 

11/8/04 

36 

64% 

79437 302 

-0.58 

1.93 

L3b 

2/17/05 

3/24/05 

35 

67% 

79 778 747 

-1.37 

1.95 

L3c 

5/20/05 

6/23/05 

34 

65% 

75 890428 

-3.44 

2.46 

L3d 

10/21/05 

11/24/05 

34 

63% 

71662 990 

-4.06 

2.56 

L3e 

2/22/06 

3/28/06 

34 

65% 

74 789 938 

-4.13 

2.10 

L3f 

5/24/06 

6/26/06 

33 

62% 

71 172 594 

-3.72 

2.27 

L3g 

10/25/06 

11/27/06 

33 

61% 

70179 600 

-4.07 

2.45 

L3h 

3/12/07 

4/14/07 

33 

63% 

73 058 093 

-3.93 

2.23 

L3i 

10/2/07 

11/5/07 

34 

59% 

68 161 111 

-3.86 

2.38 

L3j 

2/17/08 

3/21/08 

33 

63% 

71 872 084 

-3.91 

2.57 

L3k 

10/4/08 

10/19/08 

15 

59% 

29447 798 

-4.13 

2.66 

L2d 

11/25/08 

12/17/08 

22 

54% 

40 865 517 

-0.11 

5.36 

L2e 

3/9/09 

4/11/09 

33 

39% 

44 606427 

2.65 

8.79 

L2f 

9/30/09 

10/11/09 

11 

35% 

12 844 976 

1.69 

7.03 

Laser 1 avg 




73% 


1.77 

1.99 

Laser 2 avg 




51% 


1.76 

4.80 

Laser 3 avg 




63% 


-3.38 

2.32 

Average all 




60% 


-1.38 

3.13 


laser footprints from each track fall within the DEM bound- 
aries. In 2009, we resurveyed the salar de Uyuni to check 
for topographic change that might impact ICESat elevation 
validation and found that the DEM - whose error we esti- 
mate to be less than 2.3 cm root mean square (Borsa et al., 
2008) - had risen by an average of 2.5 cm (Brunt et al., 2009). 
Since we have no information about the nature of the surface 
change between 2002 and 2009 other than the small change 
in our DEMs, our best estimate of the actual surface at any 
intermediate epoch is a linear interpolation between the two 
DEMs. For the analysis used in this paper, we account for the 
effect of the surface change by linearly interpolating (node- 
by-node) between the 2002 and 2009 DEMs to the date of 
each ICESat pass over the salar, creating a slightly different 
reference DEM for each track in each campaign. 

At the salar de Uyuni, intercampaign biases for the latest 
release of the ICESat data (R633) ranged over 10 cm, with el- 
evation biases of up to 17 cm between repeated tracks within 
a single campaign (see Fricker et al., 2005 for a summary 
of our methods). These values are of similar magnitude to 
estimates by other investigators in different locations (Shu- 
man et al., 2009; Siegfried et al., 2011; Urban et al., 2012). 
The salar de Uyuni is an ideal validation site - high-elevation 
(smaller tropospheric delay correction) with negligible cloud 
cover (little or no multiple scattering) and no topography (lit- 
tle or no elevation impact from pointing errors) - so we ex- 


pected higher accuracy in our elevation recovery than we ob- 
served. At the same time we realized that with such a large 
range of observed biases, we had an opportunity to use the 
salar de Uyuni to search for candidates for the unidentified 
error sources still affecting ICESat elevations. 

2.4 Correlations between transmit pulse parameters 
and ICESat elevations 

For this study, we undertook a systematic examination of 
the elevation impact of a number of ICESat instrument pa- 
rameters, motivated by observations made by our group 
and other investigators that some of these parameters var- 
ied systematically from campaign to campaign (e.g., Fricker 
et al., 2005; Shuman et al., 2009). We hypothesized that 
at the salar de Uyuni we would be able to observe corre- 
lations between these parameters and the elevation biases 
still remaining after improvements in ICESat orbit determi- 
nation, pointing and ranging over the mission lifetime (see 
http://nsidc.org/data/icesat/past_releases.html). Although we 
recorded and tracked instrument and environmental informa- 
tion as part of our validation activities, we had not previously 
looked for quantitative correlations between these parame- 
ters and the ICESat elevation biases over Uyuni. 

Our analysis sought to explain variability in the eleva- 
tion misfits (the differences between the ICESat and DEM 
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Fig. 1. The reference DEM used in this study, located on the salar de 
Uyuni in Bolivia. Topographic relief across the 45-by-54km DEM 
is only 80 cm, making this region of the salar one of the flattest 
natural surfaces on Earth. ICESat tracks 85 (descending) and 360 
(ascending) cross the DEM and are used for range validation. 


elevations) for all 8371 valid ICESat returns over the salar 
de Uyuni DEM, without regard to their chronological or- 
der. We regressed these misfits against a number of instru- 
ment/waveform parameters and found significant non-zero 
correlation between the misfits and (1) transmit pulse skew- 
ness, (2) transmit pulse eccentricity, (3) transmit gain, (4) 
transmit pulse energy, and (5) receive pulse energy. In the 
case of transmit pulse skewness, visual examination of the 
scatterplot between skewness and misfit (Fig. 2a) shows that 
the two are linearly correlated, with stronger correlations for 
individual campaigns than for the entire data set. Quantita- 
tively, the linear Pearson correlation coefficient R between 
skewness and misfit is 0.30 (and statistically significant) for 
the entire data set, with values for individual campaigns that 
reach 0.64 for L2b and L3c. Since parameters 1-2 are re- 
lated to the transmit pulse shape, and parameters 3-5 are di- 
rectly or indirectly related to transmit pulse amplitude, we 
concluded that characteristics of the transmit pulse were af- 
fecting ICESat range determination and were able to identify 
a potential mechanism for this effect in the ICESat Range 
Algorithm Theoretical Basis Document (ATBD) (Brenner et 
al., 2003). 

ICESat’s transmit and return pulses are recorded as his- 
tograms of energy versus time, with each histogram bin span- 
ning 1 ns (equivalent to 15 cm in range). Following conven- 
tion, we refer to these histograms as laser waveforms. ICESat 
Level 1 data post-processing identifies the times associated 
with reference points on the transmit and return waveforms 
and differences the two to obtain the pulse travel time, which 
is then scaled to obtain the range from ICESat to the surface. 


Fig. 2a. Scatterplot of ICESat elevation misfits versus the transmit 
pulse skewness for each shot, with the linear correlation between 
the two indicated by the black line. The Pearson correlation coefli- 
cient R is 0.30 for the whole data set, with higher coefficients for 
most of the individual campaigns - up to 0.64 for L2b and L3c (see 
campaign color code at right of plot). 


Two types of reference point are used in ICESat processing: 
the centroid of the waveform (yielding time Cj for the trans- 
mit waveform or time Cr for the return waveform) and the 
peak position of the Gaussian fit to the waveform (times Gj 
or Gr). For consistency, travel times should be calculated us- 
ing the same type of reference point on both the transmit and 
the return waveforms (e.g., T —C^ — Cj or T — Gr — Gj). 
We noticed, however, that two different reference points were 
specified in a figure in the ICESat Range ATBD, with the 
range calculated from the Gaussian-Centroid time difference 
Gr — Ct (see Fig. 3). This is problematic because any dif- 
ference between the transmit centroid and transmit Gaussian 
peak location also appears in the return centroid and Gaus- 
sian. While this has no effect when pulse travel time is calcu- 
lated between centroids or between Gaussians, when mixed 
reference points are used it results in a range error that prop- 
agates through the ICESat geolocation process to yield an 
elevation error of opposite sign. 

2.5 An ICESat range error: the Gaussian-Centroid 
(G-C) offset 

Although Gaussian-minus-centroid (or G-C) timing was not 
intended to be used for ICESat range determination, NASA’s 
ICESat Science Computing Facility confirmed that G-C tim- 
ing was implemented through data release R633 for all 
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Fig. 3. ICESat range determination illustration, adapted from Bren- 
ner et al. (2003). The transmit and return waveforms are shown in 
blaek, and the Gaussian fits to those waveforms are shown in green. 
The eentroids of the transmit (Cx) and return (Cr) waveforms are 
indieated by the solid blaek vertieal lines, and the peak loeations of 
the Gaussian fits to the transmit (Gj) and return (Gr) waveforms 
are indieated by the dotted green vertieal lines. The ICESat range 
determination algorithm that was implemented through data release 
R633 used the time differenee between the return Gaussian peak 
and the eentroid of the transmit waveform (multiplied by the speed 
of light c and divided by 2 to get one-way range), whieh introduees 
a range error equal to (Gj - Cj). 


Fig. 2b. Same as Fig. 2a, but using G-C-eorreeted ICESat eleva- 
tions. R is now only —0.09 for the whole data set, indieating that 
the G-C eorreetion has removed most of the spurious eorrelation 
between transmit pulse skewness and ICESat elevations. 

ICESat products except GLA14 (land surface elevations). As 
discussed earlier, in the case of GLA14 the expectation was 
that over land the reflecting surface would be complex and 
the elevation of the laser footprint location would be best 
summarized by the centroid of the return waveform. GLA14 
elevations were therefore implemented using centroid-to- 
centroid (Cr — Cj) timing and are not impacted by the G-C 
offset. For GLA06, GLA12, GLA13, and GLA15, where re- 
turn waveforms were expected to be primarily single-peaked 
and Gaussian-shaped, a Gaussian fit was used to determine 
the time of the return waveform and was (incorrectly) differ- 
enced with the centroid of the transmit waveform. 

Fortunately, the range error in these four data products due 
to G-C timing (henceforth the “G-C offset”) can be exactly 
reproduced using the transmit pulse parameters available in 
the ICESat GLA05 data product: 

G-Coffset(m) = (GLA05.d_parmTr(2) 

- GLA05.d_locTr).c/2, (1) 

where d locTr is the time in ns corresponding to the transmit 
waveform centroid (parameter names are those used in the 
GLA05 data product), d_parmTr(2) is the time in ns of the 
peak of the Gaussian fit to the transmit waveform, and c is 
the speed of light defined as 0.30 mns“^. The elevation im- 
pact of the G-C offset can be removed by adding the offset 
from (1) directly to ICESat elevations, a step we will refer to 
as the “G-C correction” hereafter. Alternatively, investigators 


can apply the G-C correction from files provided by the Na- 
tional Snow and Ice Data Center (http://nsidc.org/data/icesat/ 
correction-to-product-surface-elevations.html). Technically, 
the G-C offset/correction should also be scaled by the cosine 
of the laser pointing angle measured from nadir, but since this 
angle is rarely over 3°, the scaling is < 1 mm and is negligible 
in most cases. 

After calculating and applying the G-C correction to ICE- 
Sat returns from the salar de Uyuni, we revisited our analysis 
in Sect. 2.4 and confirmed that the G-C offset was responsi- 
ble for most of the observed correlation between ICESat el- 
evations and the listed waveform/instrument parameters. As 
Fig. 2b illustrates for transmit pulse skewness, we observed 
that correlations drop to nearly zero after the G-C correction 
is applied to ICESat elevations. 

2.6 G-C offset characteristics for ICESat’s three lasers 

We took the entire ICESat data set and used Eq. (1) to calcu- 
late the G-C offset for all laser shots that resulted in a valid 
surface return. Ordering the pulses sequentially by laser shot, 
we found significant and systematic differences in the G-C 
offset over the life of each laser and between the three lasers 
(Fig. 4; Table 1). 

During its short 56-day life. Laser 1 exhibited a mean 
G-C offset of 1.77 cm, a G-C offset standard deviation of 
1.99 cm, and little change in G-C offset behavior over time 
(Fig. 4a). Laser 2 had almost the same mean offset as Laser 
1 (1.76 cm), but more than double the G-C offset standard 
deviation (4.80 cm; Fig. 4b). In fact. Laser 2 exhibited a 
four-fold increase in G-C offset standard deviation from L2c 
to L2f, which was associated with low (sub-20 mJ) and de- 
clining transmit energy. Increased variance is expected with 
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a) Laser 1 b) Laser 2 
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cumylBlIve numbef ol elwte 

Fig. 4. The G-C offset, transmit energy, and transmit gain for Laser 
1 (a), Laser 2 (b), and Laser 3 (c), ordered sequentially by the eu- 
mulative number of shots for eaeh laser. The red line in the top plot 
is the 10 000-shot moving average of the G-C offset. ICESat eam- 
paigns are named in blue at the top of the plots and are delineated 
by blue lines marking the end of eaeh eampaign period. The seatter 
in the G-C offset grows as transmit energy drops, espeeially below 
20 mJ. 


low transmit energy because the accompanying decrease in 
the signal-to-noise ratio of the laser waveforms (once the 
transmit gain can no longer be increased to compensate) de- 
grades the precision of centroid determination and Gaussian 
fitting (Pricker et al., 2005). In addition, there were signif- 
icant changes in the moving average of the G-C offset (red 
line in Fig. 4), including (1) a 5 cm drop at the end of L2a 
that is associated with an instantaneous 10 mJ fall in transmit 


ICESat 6-C offset, campaign means and trends 



Fig. 5. G-C offset campaign means and trend estimates for selected 
data spans from Table 2, Trend A. The box symbols indicate the 
mean values of the G-C offset for each campaign, the error bars 
show the offset standard deviations, and the different lines are G- 
C offset trends estimated for different data periods using inverse- 
variance weighing. The G-C offset has a bigger impact on elevation 
trends when campaigns at the end of the mission are excluded. 

power and is probably due to the failure of one of the diode 
pump bars, (2) a large negative spike at the beginning of L2b 
during the period ICESat was in sun acquisition mode, and 
(3) an inverse correlation with transmit power beginning in 
L2d, which was when transmit energy fell and remained be- 
low lOmJ. 

By contrast. Laser 3 showed much more stable G-C off- 
set behavior, due in part to its higher transmit energies 
(Fig. 4c). Overall, Laser 3 had a lower mean G-C offset value 
(—3.38 cm) than the other two lasers and a relatively constant 
standard deviation (2.32 cm) whose value was close to that of 
Laser 1 . The G-C offset moving average started around 0 cm, 
jumped 3 cm in the middle of L3a (when the laser tempera- 
ture was increased from 13.8 to 16.0 °C), and then gradually 
dropped by 4 cm through L3d and remained consistent up to 
the end of L3k. 

3 Results and discussion 

3.1 Trend in the G-C offset over the ICESat mission 

There are two direct effects of the G-C offset: (1) it increases 
the shot-to-shot variability of ICESat elevations (especially 
at low transmit energies) and (2) it shifts the mean elevation 
for most campaigns. Of greater relevance to ice-sheet mass 
balance studies is a secondary effect due to the fact that when 
ordered in time, the changes in the mean G-C offset between 
campaigns exhibit time-correlated behavior over the mission 
period (Fig. 5) that could potentially be interpreted as real 
surface elevation change. In particular, the high values of the 
G-C offset for all three campaigns at the end of the mission 
(L2d-L2f) generate erroneously low elevations that could be 
interpreted as suddenly higher ice sheet mass-loss rates at the 
end of the mission period. 

Most published studies using ICESat data for ice sheet 
mass balance have derived surface change estimates by fit- 
ting linear trends in time through all elevations over a target 
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Table 2. G-C offset trend estimates (and formal 1 -sigma errors) using different data periods and different least-squares weighting applied 
to eaeh eampaign. Trend A uses inverse-varianee weighting derived from the G-C offset statisties in Table 1, Trend B uses uniform weight- 
ing, and Trend C uses weights derived from global ICESat sampling by eampaign. These trend estimates represent the impaet of the G-C 
eorreetion on d/z/dt, indieating that all ICESat-derived iee elevation ehange estimates should beeome more negative onee the G-C offset is 
removed. 


Begin 

campaign 

End 

campaign 

G-C Trend A 
(cm yr“^) 

G-C Trend B 
(cm yr“ ^ ) 

G-C Trend C 
(cm yr“ ^ ) 

L2a 

L3i 

-2.13 

±0.53 

-1.90 

±0.48 

-1.90 

±0.51 

L2a 

L3j 

-1.88 

±0.47 

-1.62 

±0.41 

-1.76 

±0.47 

L2a 

L3k 

-1.64 

±0.41 

-1.39 

±0.36 

-1.55 

±0.42 

L2a 

L2d 

-1.54 

±0.40 

-1.00 

±0.32 

-1.17 

±0.38 

L2a 

L2e 

-1.48 

±0.40 

-0.62 

±0.29 

-1.05 

±0.37 

L2a 

L2f 

-1.38 

±0.40 

-0.36 

±0.26 

-0.92 

±0.36 


area (Schenk and Csatho, 2012 offer an alternative method- 
ology that accounts for non-linear features in elevation time 
series). Since we wanted to understand how the G-C offset 
might have impacted these mass balance estimates, we es- 
timated the trend of the G-C offset trend over several time 
periods. First, we performed a linear regression against time 
of the global mean G-C offsets for campaigns L2a to L2f 
(the entire mission), with inverse-variance weights for each 
campaign calculated from the G-C offset standard deviations 
(all data are from Table 1). This choice of weighting low- 
ers the contribution of high- variance data points, which in 
this case are the three Laser 2 campaigns at the end of the 
mission. For the L2a-L2f period, we obtained a trend of 
— 1 .4 zb 0.4 cm yr“ ^ , which is statistically different from zero 
at the 3-sigma level. Since many published studies do not in- 
clude the later ICESat campaigns in their analysis, we also 
report trends for other data periods (Table 2, Trend A). As 
Table 2 and Fig. 5 show, the fewer campaigns used at the end 
of the ICESat mission, the more negative the trend of the G-C 
offset: up to —2.1 ± 0.5 cmyr“^ in the case of data spanning 
only the period between L2a to L3i. 

3.2 Potential impact of the G-C offset on ICESat 
elevation trends 

The negative trend in the G-C offset contributes an erroneous 
positive trend in ICESat elevations that could be interpreted 
as real surface change. However, the trends we calculated 
in Sect. 3.1 are not our best estimate of the actual effect of 
the G-C offset on ICESat-derived elevation change because 
of differences in linear regression weighting, as we discuss 
below. 

ICESat investigators typically estimate volume change by 
integrating many independent elevation trend estimates over 
a region of interest. If ICESat elevation data are not corrected 
for the G-C offset, the G-C offsets on all the individual laser 
shots will propagate through the linear regression used to es- 
timate elevation trends and will introduce an error. Impor- 
tantly, this trend error is highly dependent on the weighting 


applied to the elevations in the regression (see Appendix A). 
Since most investigators use constant weighting for all re- 
turns (e.g.. Shepherd et al., 2012), we repeated the linear 
regressions from Sect. 3.1 using a uniform weight on each 
campaign of 1/4 cm“^ (corresponding to a 2.0 cm standard 
deviation). The result was less negative G-C offset trends 
(Table 2, Trend B), with the biggest change for the longest 
data periods (e.g., for L2a-L2e and L2a-L2f). Using a differ- 
ent value for the uniform weight will not change these trend 
estimates, although it will change the formal error. 

In addition, the number of ICESat measurements in each 
campaign (i.e., spatio-temporal sampling) also affects how 
the G-C offset will affect elevation trend estimates. For ex- 
ample, campaign L2f had fewer returns than any other cam- 
paign because it was only 1 1 days long and because its low 
laser transmit energy resulted in a low percentage of valid re- 
turns being recorded due to attenuation by clouds. There are 
relatively few L2f elevations included in the many thousands 
of elevation trend estimates made over an ice sheet, which 
means that campaign L2f elevations (and therefore the G-C 
offset for L2f) are typically underweighted in the overall ice 
sheet elevation trend. 

We estimated the impact of sampling on the G-C trend by 
modifying the uniform- variance regression weights above to 
include a term for the relative number of returns expected 
for each campaign (Appendix A). The resulting trend esti- 
mates (Table 2, Trend C) are more negative than they would 
be if sampling were ignored. This is due to the lower weight- 
ing of the L2d-L2f campaigns (whose low energies resulted 
in fewer returns than in earlier campaigns), which reduced 
the impact of their relatively high G-C offsets on the G-C 
trend estimate. These trends are our best approximation of 
the impact of removing the G-C offset from ICESat eleva- 
tions, although differences in relative ICESat sampling from 
the global campaign averages (due to regional effects or data 
editing protocols) are likely to cause the impact of the G-C 
offset in specific cases to vary from our estimates. 
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3.3 Power spectrum of the G-C offset 

To determine if there is regional variation in the G-C off- 
set, we examined the G-C offset power spectrum to look for 
temporal correlations that could map into spatial patterns. 
For this analysis we could not use G-C offsets calculated 
from the centroid and Gaussian parameters in the GLA05 
data product, since shots without a valid return (see Table 1) 
did not undergo Gaussian fitting during ICESat data process- 
ing. Spectral estimation requires continuous time series, so 
we retrieved transmit waveform records from the GLAOl 
data product, estimated Gaussian fits for every ICESat trans- 
mit pulse, and recalculated the G-C offset using Eq. (1) 
and the original GLA05 centroids. The G-C offsets calcu- 
lated this way deviate from the GL AOS -derived offsets by 
0.0 zb 2.5 mm (1-sigma) overall and are continuous over each 
campaign. 

We made G-C offset power spectral density (PSD) esti- 
mates for each of the 18 ICESat campaigns using a single- 
window Fast Fourier Transform (FFT) with a Hamming ta- 
per. In all cases, the power spectra are flat at frequencies 
higher than 10“^ Hz, except for narrow spectral peaks at 
0.000173 Hz (the orbital frequency of ICESat) and various 
multiples thereof (Fig. 6). If there were any jumps or ramps 
in the G-C offset time series during a campaign, these are 
manifested as a ramp in the spectrum at frequencies lower 
than 10“^ Hz (e.g.. Fig. 6, left). What these results indicate 
is that the G-C offset generally behaves like white noise over 
periods shorter than a day (about 15 full orbits), with super- 
imposed sinusoidal variability over distances < 1 orbit. Al- 
though this implies some degree of along-track correlation, 
because successive ICESat ground tracks All in between ear- 
lier tracks, changes in the G-C offset are almost uniformly 
distributed over a broad area such as an ice sheet. The global 
characteristics of the G-C offset should thus be a good first 
approximation to how it manifests in any regional analysis. 

3.4 Impact of the G-C Offset on ICESat d/^/dt estimates 
over ice shelves 

In order to test our assumptions about the impact of the G-C 
offset on ICESat trend estimates, we did an analysis of data 
over the Ross and Filchner-Ronne ice shelves in Antarctica. 
These ice shelves are large, relatively featureless, and are a 
good analogue for the low-relief interiors of the ice sheets. 
For each ice shelf we estimated elevation trends for cam- 
paigns L2a-L2f using ICESat data with and without the G-C 
offset. For this analysis, we followed the standard approach 
of simultaneously estimating planar slopes (d/i/dv, dhidy) 
and temporal trends (dhldt) for nearby footprints along seg- 
ments of ICESat repeat tracks using least-squares estima- 
tion with unit weighting (e.g.. Smith et al., 2009; Gardner et 
al., 2013). We used saturation-corrected and tide-corrected 
elevations from the ICESat GLA12 data product, “retided” 
the elevations using a more accurate tide-model (Padman et 


al., 2002; Fricker and Padman, 2006), and estimated dh/dt 
only for tracks that were repeated during four or more cam- 
paigns. Finally, for each ice shelf and data set, we calculated 
a latitude-weighted average dhidt value by applying an em- 
pirical function of the satellite orbit convergence towards the 
geographical poles (Gardner et al., 2013). This provides a 
weighted estimate of dhidt that is consistent with gridding 
methods, although without the need for spatial interpolation. 

For both ice shelves, the impact of the G-C correction 
was that the average dhidt value became more negative, on 
the order of 0.6cmyr“^ For the Ross Ice Shelf, the aver- 
age dhidt value changed by — 0.68cmyr“^ after applying 
the G-C correction (from -h0.16 cmyr“^ to — 0.52 cmyr“^). 
For the Filchner-Ronne Ice Shelf, the average dhidt value 
changed by —0.49 cmyr“^ after applying the G-C correction 
(from -hi. 78 cmyr“^ to -hi. 29 cmyr“^). The reason the im- 
pact of the G-C correction is different on each ice shelf is 
because the spatio-temporal sampling in each campaign is 
different. Most relevant is the undersampling of the Filchner- 
Ronne Ice Shelf in campaign L2a relative to the Ross Ice 
Shelf (see Table 3), which slightly flattens the G-C correc- 
tion trend. This can be understood by observing that reduc- 
ing the weight of L2a in Fig. 5 will mitigate the impact of its 
high G-C offset and thus flatten the slope of the linear fit to 
the G-C offsets. Sampling in campaigns L3g and L3i is also 
different on the two ice shelves, but these differences matter 
less because L3g and L3i are close to the center of the time 
series and because their mean G-C offset values are similar 
to those of nearby campaigns. 

The G-C correction significantly changed the mass bal- 
ance of the two ice shelves we examined. Because of hy- 
drostatic equilibrium, a given change in surface elevation 
equates to about nine times more change in ice thickness, 
greatly magnifying the mass balance impact of systematic 
measurement errors. In the case of the Ross Ice Shelf, the 
— 0.68cmyr“^ change in average dhidt from the G-C cor- 
rection implied a mass balance correction of —29 Gt yr“^ af- 
ter accounting for hydrostatic equilibrium. The mass balance 
correction for the Filchner-Ronne Ice Shelf was 18 Gtyr“^ 

For comparison, we also approximated the impact of the 
G-C correction on the ice shelf elevation trend using linear 
fits to the mean G-C offsets of each campaign. For each ice 
shelf, we linearly regressed the mean G-C offsets from Ta- 
ble 1 against time, using the campaign sampling from Table 3 
to derive appropriate weights via Eq. (A 10) (Appendix A). 
Our estimates for the impact of the G-C correction using 
this method were —0.61 cmyr“^ for the Ross Ice Shelf and 
— 0.49cmyr“^ for the Filchner-Ronne Ice Shelf, which are 
both within 1 mmyr“^ of the actual dhidt values calculated 
by applying the G-C correction and reprocessing the data set. 
This suggests that it may not be necessary to fully repeat ear- 
lier ICESat studies to determine the impact of the G-C cor- 
rection as long as the local sampling and explicit weighting 
scheme for each campaign is known. 
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Table 3. ICESat sampling of the Ross and Filchner-Ronne ice 
shelves for the elevation trend estimates in Sect. 3.4. We list the 
number of shots used by campaign over the entire ice sheet (5/ in 
Eq. AlO), as well as the relative campaign sampling weights (ex- 
pressed as a percentage of the maximum number of shots on that 
ice shelf for any campaign). There are large differences in relative 
campaign sampling between the two ice shelves for L2a, L3g and 
L3j, which results in different estimates for the elevation trend im- 
pact of the G-C offset. 


Campaign 

Ross 

# Shots used 

% Max 

Filchner-Ronne 
# Shots used % Max 

L2a 

395218 

96 

179 955 

68 

L2b 

285 648 

69 

198316 

75 

L2c 

378 640 

92 

244 073 

92 

L3a 

377 854 

91 

245 864 

93 

L3b 

335 248 

81 

202 618 

77 

L3c 

363 814 

88 

264 081 

100 

L3d 

358 902 

87 

197 763 

75 

L3e 

291486 

71 

151919 

58 

L3f 

368413 

89 

234 540 

89 

L3g 

413 425 

100 

169199 

64 

L3h 

329 989 

80 

227 046 

86 

L3i 

384517 

93 

235 922 

89 

L3j 

262 832 

64 

220 324 

83 

L3k 

160299 

39 

127 620 

48 

L2d 

189205 

46 

123 881 

47 

L2e 

261577 

63 

190 585 

72 

L2f 

104189 

25 

49165 

19 


3.5 Impact of the G-C offset on ICESat intercampaign 
biases 

There have been many independent estimates made of ICE- 
Sat intercampaign biases (e.g., Siegfried et al., 2011; Zwally 
et al., 2011; Shepherd et al., 2012; Ewert et al., 2012; Gunter 
et al., 2013; Hofton et al., 2013), many of which were di- 
rectly compared in Urban et al. (2012). There is little con- 
sensus between estimates of individual biases or bias trends: 
the seven estimates considered in Urban et al. (2012) yielded 
biases that differed by up to 20 cm for any single campaign 
and bias trends that ranged from —0.3 to +2.2 cmyr“^ over 
the L2a-L2f period. Given the different surface types, data 
locations, spatial sampling, and methodologies used in these 
estimates, it is not surprising that the estimated biases dif- 
fer. However, since the intercampaign bias trends (or their 
underlying biases) are supposed to be applied as corrections 
to all ICESat data, it is important to understand the reasons 
for and implications of the variability between the different 
estimates. 

We identify three components of the measured intercam- 
paign biases: (1) the contribution of the G-C offset; (2) bias 
due to all other instrument error sources; and (3) bias esti- 
mation errors due to actual elevation/range changes from un- 
modeled physical processes such as surface change and at- 
mospheric scattering over the selected calibration surfaces. 


The contributions of all three components will vary accord- 
ing to the data and methodology used in the bias estimation. 
For instance, we demonstrated in Sect. 3.4 how the impact 
of the G-C offset differs by a few mm yr“ ^ between two ice 
shelves on the same continent because of different spatio- 
temporal data sampling. We might also expect, for example, 
differences between land-, ocean- and ice-determined biases 
due to systematic differences in cloud cover or surface ge- 
ometry (e.g., sea-state effects). Larger differences are likely 
if biases are being estimated over surfaces that are assumed 
to be stable but are not, or over surfaces whose time variance 
is imperfectly known. 

Applying the G-C correction removes the first component 
of the intercampaign biases and thus will alter existing in- 
tercampaign bias estimates. Specifically, we can expect the 
G-C correction to decrease intercampaign bias trends by 0.92 
to 1.90cmyr“^ depending on the campaigns used (Table 2, 
Trend C). This means that the G-C offset is a significant 
contributor to the intercampaign biases, although the large 
spread of the intercampaign bias estimates suggests that sig- 
nificant residual variability will remain in all estimates af- 
ter the G-C correction - at the level of individual campaign 
biases and/or in the overall trends themselves. Validation 
using the two ICESat tracks that cross the salar de Uyuni 
DEM (Fig. 1) confirms that the impact on the fitted eleva- 
tion trend from applying the G-C correction to individual 
ICESat shots is almost the same as what we predict from 
the global analysis in Sect. 3.2 (—0.92 cmyr“^ predicted vs. 
— 1 . 1 7 cm yr“ ^ actual for L2a-L2f). However, the magnitude 
of the Uyuni intercampaign bias trend (which changed from 

0.67 ± 0.47 cmyr“^ to —0.50 ± 0.36cmyr“^) did not appre- 
ciably drop, which illustrates our point that post-correction 
errors are likely to remain despite improvements in accuracy 
and precision. 

Of importance to previous studies that used data contain- 
ing the G-C offset, recent studies have suggested (e.g., Rig- 
not et al., 2013) that applying a set of intercampaign biases 
implicitly corrects for the impact of the G-C offset, at least 
at the level of the mean elevation for each campaign. While 
technically true, we are concerned that this is not a satisfac- 
tory way to make the G-C correction. Applying empirical 
intercampaign biases may partially correct the effect of the 
G-C offset, but the variability of intercampaign biases from 
different sources suggests that this approach can introduce 
additional errors in dhidt estimates. We recommend instead 
that investigators explicitly correct for the known G-C offset 
using one of the methods described in Sect. 2.4. 

Finally, we consider whether intercampaign biases should 
be applied after making the G-C correction. Investigators 
should be aware that: 

1 . There is still no authorized or consensus set of inter- 
campaign biases. 

2. Differences between biases from different sources can 
be statistically significant. 
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G-C offset power spectrum, L3A G-C offset power spectrum. L3C 



Fig. 6. Power spectrum of G-C offset for campaigns L3a (left) and L3c (right), showing a flat spectrum beyond 10“^ Hz, narrow spectral 
peaks at harmonics of the orbital frequency of ICESat, and different behavior at lower frequencies. The main difference between the two 
campaigns is that the moving average of the G-C offset jumps abruptly in L3a while it is flat in L3c (see Fig. 4b). 


3. Bias estimates must be G-C-corrected before being ap- 
plied to G-C-corrected elevations. 

4. The biases for each campaign should be applied to in- 
dividual elevations before dh/dt is estimated. Those 
who wish to apply intercampaign bias trends directly 
to dh/dt estimates should be aware that they must cor- 
rect bias trends for sampling effects and the weights 
used in the dh/dt calculation, as discussed earlier. 

4 Conclusions and outlook 

We have reported on a range correction to the ICESat Level 
1 data that removes the effect of an erroneous travel time 
calculation in the Level 1 data (the G-C offset). The impact 
of the G-C offset on ice-sheet elevation trends (dh/dt) can 
vary substantially depending on the time span of investiga- 
tion and the data sampling in each observation campaign. 
If those factors are carefully considered, we have shown 
that it is possible to reproduce the effect of the G-C off- 
set to within a few mmyr“^ at a regional scale. Ideally, 
users should apply the G-C offset corrections to the data 
themselves based on the GLA05 data product or use cor- 
rection files provided by NSIDC (http://nsidc.org/data/icesat/ 
correction-to-product-surface-elevations.html) until a final 
release of the ICESat data becomes available. 

Additional work is still needed to characterize the eleva- 
tion errors that will remain after the G-C offset is removed. 
This work includes the revision and ultimate reconciliation 
of various estimates of ICESat intercampaign biases, which 


will change as a result of the G-C offset correction. The large 
variation in existing estimates of intercampaign biases sug- 
gests that the problem of estimating empirical errors is not 
necessarily any easier than uncovering the root sources of 
those errors. We also wonder if the formalization of the in- 
tercampaign bias as a description for otherwise unmodeled 
and persistent errors in ICESat elevations may have diverted 
attention away from the need for a more systematic and me- 
thodical effort to identify outstanding error sources in the 
ICESat data. In particular, we are concerned that a single 
“universal” set of intercampaign biases (were one to become 
available) would not be equally relevant across a range of 
studies using different spatial subsets of data or different data 
sampling/editing. 

The discovery of the G-C offset was the result of having 
access to an accurately surveyed stable reference surface that 
allowed us to unambiguously link a component of ICESat- 
measured elevation changes to an ICESat timing error. Our 
approach of correlating elevation misfits with various instru- 
ment parameters can and should be applied to a broader sam- 
ple of reference surfaces to provide a greater range of param- 
eter values, thereby increasing the diagnostic power of the 
correlation analysis. There are many active and potential ref- 
erence sites around the globe, and we believe that the satel- 
lite altimetry community should attempt to link these into a 
single virtual surface for altimeter calibration and validation. 
Future missions might consider increasing both the resources 
and expectations for these efforts. 
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Appendix A 


Linear regression of mean G-C offsets to estimate the 
G-C impact on ICESat elevation trends 


We take the chi-square merit function for the linear fit to N 
pairs of time/elevation data to be 


x\a,b) 


N-l r 


E 


(hi -h Ahi) — (a bti) 


n2 


(Al) 


where the tt are points in time corresponding to different 
ICESat campaigns, hi are the true surface elevations, Ahi 
are the corresponding G-C offsets, and <Ji are the standard 
deviations assumed for (hi -h Ahi). The maximum likelihood 
estimate of the slope (or trend) b of the linear fit is the closed- 
form solution obtained from minimizing Al (Press et al., 
2007) 
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with separate terms for the true surface volume change and 
volume change from the G-C offset. If we make the simplify- 
ing assumption that the area Aj of each evaluation cell is the 
same (most studies interpolate volume/mass estimates to the 
nodes of a regular grid) and take M to be the number of cells, 
the G-C volume change term from Eq. (A5) can be written 
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Substituting into Eq. (A6) the linear regression solution for 
bo-c from Eq. (A3), taking the values for ti and a/ to be the 
same everywhere for a given campaign (although different 
between campaigns), and distributing the summation over j 
gives 
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Since A2 is linear in h-\- Ahi, it can be rewritten as the 
sum of separate solutions for h and Ahi 
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or 

b = ^surface +^G-C- 


(A3) 


(A4) 


The contribution of the G-C offset to the trend in the data 
can therefore be considered separately from the contribution 
of the surface elevations themselves, although the weighting 
will be the same in both cases. 

To scale these conclusions from a single evaluation cell to 
an entire study area, we note that the ice volume change for 
a given region is calculated by summing the contribution of 
many independent trend estimates 
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where 
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is the mean G-C offset value for the campaign indexed by /, 
and 


EA 


XI A 




X A A/i/ 




(A9) 


is the maximum likelihood estimate of the G-C offset trend 
from the mean G-C offsets Ahi. 

Estimates of Ahi are given in Table 1 (mean G-C offset by 
campaign), which our analysis in Sect. 3.3 suggests should 
be valid anywhere on Earth. Equation (A9) is evaluated by 
estimating the linear trend Ahi ordered in time (using what- 
ever Oi were chosen for the surface change analysis) and the 
result is used in Al to estimate the impact of the G-C offset 
on ICESat elevation change trends assuming identical sam- 
pling of the various ICESat campaigns. 

Finally, we consider the impact of sampling, whereby the 
number of valid ICESat returns can be different for differ- 
ent campaigns. In Eq. (A3), this is manifested by having 
some campaigns with a large number of elevations/offsets 
and some with few (or even no) elevations/offsets, which 
implicitly down-weights campaigns with fewer samples. We 
would like to have a way of representing sampling when we 
fit linear trends to the G-C offset means in Eq. (A9). While 
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it is beyond the scope of this paper to offer a proof, boot- 
strap analysis will confirm that a weighting that accounts for 
non-uniform sampling is 


±TA 

Si N ’ 


(AlO) 


where the St are the number of samples for a given campaign 
in the study area and N is the total number of campaigns (as 
in Eq. Al). This formulation roughly preserves the value of 
the formal error on the slope estimate obtained from uniform- 
variance weighting. In the uniform-variance scenario typical 
for ICESat, what matters for the trend estimate is the rela- 
tive number of samples between campaigns, not the absolute 
numbers. 
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